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We propose a model based on coupled multiplicative stochastic processes to understand the dy- 
namics of competing species in an ecosystem. This process can be conveniently described by a 
Fokker-Planck equation. We provide an analytical expression for the marginalized stationary dis- 
tribution. Our solution is found in excellent agreement with numerical simulations and compares 
rather well with observational data from tropical forests. 
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I. INTRODUCTION 

One of the the most widespread quantities employed in 
Ecology to describe the biodiversity in a given ecosystem 
is the distribution of species abundance. In operational 
terms it can be defined as the histogram of the number of 
species (in a well defined temporal and geographical con- 
text) consisting of a generic number of individuals, or, 
from a more theoretical perspective, as the probability 
that a generic species is composed by a certain number 
of individuals. Data collected in different locations sug- 
gests that the relative species abundance distributions 
show a certain degree of similarity Q. To elucidate the 
causes that determine the shapes of these distributions 
and therefore their similarity is a problem of the utter- 
most importance and not only of theoretical nature: to 
understand the motives that influences the relative rar- 
ity or commonness of different species can be of great 
help in determining policies for the conservation of the 
endangered ones. 

The first studies on this subject can be dated back to 
the '40 and are due to Fisher |2] and Preston Their 
works were focused on finding distributions that could fit 
well particular data set in an empirical way. In particular, 
Preston |3( argued that the probability of finding species 
with a certain number of of individuals x should be log- 
normal distributed, while Fisher proposed a function 
of the form e~ ax fx, with a « 1, the so-called Fisher log 
series. 

Later, MacArthur pj firstly pointed out that similar 
distributions are found in very different ecosystems, sug- 
gesting that the shape of such distributions is to a large 
extent determined by very basic, general and ecosystem- 
independent mechanisms. This in turn hinted to the pos- 
sibility to predict the shape of such distributions with 
simple and general models, without taking into account 
too many specific details of the ecosystem under consid- 
eration. Several models have been proposed that spoused 
this view 0, IE 111- Many of them restrict to model- 
ing a single ecological community, a collection of simi- 



lar species that feed on the same pool of resources in a 
local area. This definition implies that species belong- 
ing to the same community interact mainly in a com- 
petitive way: in particular, there are no prey-predator 
relationships among them. The particular case of single 
ecological community can be framed in the wider context 
of a neutrality hypothesis. The concept of neutrality was 
firstly introduced in the framework of a biomolecular evo- 
lution theory by Kimura 7] , and then extended to other 
fields of biology. In the words of Hubbell 0, an eco- 
logical theory can be considered neutral when "...treats 
organisms in the community as essentially identical in 
their per capita probabilities of giving birth, dying, mi- 
grating and speciating. This neutrality is defined at the 
individual level, not the species level ..." 

The question whether there exist ecological communi- 
ties satisfying this assumption is still rather controversial 
@ , therefore it is crucial to understand what are the con- 
sequences of this zero order hypothesis 9] . In the context 
of a neutral hypothesis it is reasonable to describe the 
number of offspring to which any given individual gives 
place to as a stochastic variable. As a consequence, the 
number of individuals in a species at a given time can be 
regarded as a multiplicative random process. 

Here we present a model aimed at reproducing the fea- 
tures of species abundance distributions under a minimal 
set of assumptions: neutrality and the possibility to de- 
scribe the birth process as a multiplicative random pro- 
cesses. The model translates in a Fokker-Plank equation 
for the species abundance distribution and is amenable 
to an analytical treatment. The solutions found are com- 
pared with the experimental data we avail of. The shape 
of these solutions depend on one parameter, and give in 
the two limiting cases both a lognormal-like curve and 
the Fisher log series. The paper is organized as follows. 
In the second section we will present the model and com- 
ment the assumptions made. In the third one we will 
take the continuum time limit of our model, and will 
provide an analytical solution for the marginalized sta- 
tionary probability distribution function. In the last two 
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sections we compare our results with the experimental 
data and comment them. 



II. THE MODEL 



The number of individual of each species will be typi- 
cally large, so we apply the central limit theorem to the 
sum of random variables in this equation, obtaining the 
following model: 



Let us consider an ecological community consisting of a 
fixed number, s, of species. According to MacArthur and 
Wilson theory of island biogeography 0] j the number 
of species in a community approaches a dynamical equi- 
librium between immigration, speciation and extinction. 
We assume that we can neglect the fluctuations around 
this equilibrium value: in our model, when a species go 
extinct, it is immediately replaced by another one. We 
also assume that the net effect of the competitive inter- 
action between species in the community is just to keep 
also the total number of individuals in the community 
fixed: the resources available are enough to support just 
TV individuals across all the species. This last assump- 
tion implies that the populations of the species undergo 
a zero sum dynamics. This hypothesis is well confirmed 
by experimental data 0, E3; at the end of section III 
we will show that relaxing these constraints does lead to 
similar conclusions in the large N limit. We introduce 
the s variables x\ , representing the population of the «-th 
specie at (discrete) time t, with the condition: 



= N 
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where A and a are the mean value and the r.m.s.d. of the 
distribution P(X), and the £'s are uncorrelated gaussian 
variables with zero mean and unit variance. 

It is worth noting the relation between our model and 
the multiplicative process introduced by Kesten in [Tlj . 
Kesten studied random multiplicative processes of the 
form X t +i = XtX t + b t , where X t is the variable and both 
A and 6 are random variables. He found that, depending 
on the mean value of A and on the boundary conditions, 
one retrieves a lognormal or a power-law regime. Models 
for ecology and economics based on this kind of processes 
were proposed by Sornette |12| and Solomon |13|| . In our 
model the number of individuals of different species can 
be thought as following coupled Kesten-like processes. 
The coupling is a consequence of the constrain that keeps 
fixed to ./V the number of individuals in the community 
and that is enforced in equation (1) by the factor N and 
by the denominator. 
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Let -P(A) be the probability that an individual in the 
community has A offspring during one time step. Here 
neutrality plays a key role: our assumption implies that 
P(X) is the same for all individuals. The population of 
the i-th species evolves according to the following equa- 
tion: 
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where [ ] means the integer part. We are assuming that 
the existence of species with a non integer number of 
individuals is not too drastic. This might lead to round- 
off problems only for rare species. At each time step 
(generation) we just sum the number of offspring of ev- 
ery individual belonging to that species, and then add a 
small quantity b. This quantity becomes relevant only 
for small x%, and this describes the behavior of species 
near their extinction threshold. We are assuming that 
the net effect of extinctions, immigration and speciation 
can be modeled in a simple way with this term, whose 
effect is to force the Xi's to be greater than zero. Indeed, 
for 6 = 0, our system admits an absorbing state with 
only one Xi equal to N and the others equal to 0, the so- 
called monodominance Q. Notice that species are only 
coupled through the denominator, that simply preserves 
the normalization condition. 



III. THE CONTINUUM LIMIT 

In order to obtain some analytical result, we do the 
continuous time limit of this model, by introducing the 
time interval dt in the following way: 



A 
6 
a 



1 + Xdt 
b dt 
a dt 



(3) 



By means of this substitution, our model becomes: 



„t+dt 



x\ + dt(\x\ + (Ty/X\Q + b) 



(4) 



Expanding the denominator and using the fact that 
J2j x j — N, we get the Langevin equation: 



where: 



Xi = fi(x) +^2Bi j (x)£ j 

3=1 



fi(xi) = 6(1 - —Xi) 
Bij(x) = [Sij - jf)V^J 



(5) 



(6) 
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The Fokker-Planck equation associated to this 
Langevin equation is : 



with D = \ and: 



-fiP(x, t)+D^2 d 3 (gMx)P(x, t)) 

3 

(7) 



9ij 



(x) = gji(x) = BikBjk = (dij - -^) 



(8) 



We search for a solution of this equation satisfying de- 
tailed balance (i.e. P st fi — D'Yljdj(gijP st )). Defining 
the marginalized probability distribution function: 



p(x) 



Y\d Xj P st {x) 



(9) 



we can easily obtain an equation for p{x) 

d 
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This equation can be easily solved, giving: 



p(x) oc x 



0-1 



x\P( s —k)- 1 
iV/ 







D 
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Notice that this distribution correctly shows the mon- 
odominance behavior 6(0) or 6(N) in the limit j3 — * 0. 



Finally, if we fix [i 
obtain: 



p(x) 



JV 



, in the limit for N — > oo we 



(12) 
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In fig. 1 we plot simulation of the stationary p.d.f. for 
various value of the parameter /3, and check the validity 
of O- 

Instead of having a system of stochastic differential 
equation, it is possible to take into account the interac- 
tion of a species with the ecosystem in an averaged way. 
Let us consider the Langevin equation: 



b(t) = b + Xx-^fx + Dy/x£ 



(13) 



where the parameter 7 takes into account the effect of 
competition. In order to have normalizable solutions, we 
have to require that 7 > A. When this condition holds, 
it is straightforward to show that the stationary p.d.f. 
satisfying detailed balance is the same as (I12f> . with /1 = 
— (A— j)/D. Notice that in this case, the detailed balance 
solution is exact; it is also remarkable that the stationary 
distribution (|12|1 can be achieved without fixing neither 
the number of species, nor the number of individuals. 




FIG. 1: simulation of marginalized stationary p.d.f. for vari- 
ous values of the diffusion coefficient D, compared with the- 
oretical curves. For all curves 6=1, A = ls = 100, iV = 10 9 
Curves are binned linearly with binning size 6x = 10 4 . Notice 
that as D increases the curve approaches the Fisher log series. 



IV. COMPARISON WITH EXPERIMENTAL 
DATA 

Among the most reliable data on single-trophic species 
distribution of species abundance are tropical forest cen- 
sus |15). In order to make a coarse graining, a Preston 
plot is used: data are collected via a logarithmic binning 
in base 2, and species at the edge between two consecu- 
tive binning are equally divided between them. Since we 
have a continuous probability density, we compared the 
histogram with the integral over the bins of the distribu- 
tion with the experimental data, and made a least-square 
fit of the parameters (5 and /i, plus the normalization. We 
found a good agreement of our predicted curve with the 
histogram; in FIG. 2 it is shown the comparison between 
our solution and the lognormal. Notice that the two dis- 
tributions have the same number of fitted parameter. It 
would be interesting to compare our distribution with 
data collected form other kind of ecosystems, and to try 
to clarify the dependence of our free parameter (3 from 
ecological quantities like the immigration pressure, the 
speciation rate and the extinction threshold. 



V. DISCUSSION AND PERSPECTIVES 

The model we introduce admits a family of stationary 
p.d.f. depending on the parameter f3. This parameter 
fully determines the shape of the distribution: for << 
1 one recovers the Fisher log series, while for (3 large, 
one obtains a lognormal-likc distribution. As we already 
pointed out, both these distributions are well known in 
the population biology literature as possible candidate to 
be the 'right' distributions found in nature. 

There is some analogy between our model and the 
Kesten process. Indeed, also the Kesten process ad- 
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— Pasoh Preston plot 
least-square fit 

- * Lognormal fit 




mits two different regimes, one lognormal and one with 
a power law tail. The main differences is that in our 
case the multiplicative random process is applied to the 
square root of the variables, rather than to the variable 
itself. As a consequence, in the Kesten case, the expo- 
nent of the power law tail of the stationary distribution is 
always greater than one, while the small (3 regime of our 
system is characterized by a power law tail over many 
decades, with an exponent that is always less than 1: the 
cutoff due to the conserved number of individuals ensures 
the normalization of these long-tailed distributions. 

It is remarkable that our distribution is the same found 
in studies made by Kerner in the '50 | 16| on the invari- 
ant measure in a system of Lotka-Volterra equations with 
purely asymmetric couplings. In that works the interac- 
tions are only of predator-prey type, and the system is 
deterministic, while we are considering a stochastic sys- 
tem with purely competitive coupling. The discover of 
the same distribution in such different models suggests 
that it might exist some deeper and more general mecha- 
nism determining the statistical behavior of ecosystems, 
regardless of the type of interactions among species. 



FIG. 2: fit of Barro Colorado Island and Pasoh specie abun- 
dance data - Preston plot Q. Comparison between our so- 
lution and lognormal. Fitted value of the parameters of our 
distribution are /3 = 0.23 and ^ = 0.010 for the BCI; /3 = .37 
and \x = 0.015 for Pasoh. In absence of an objective estimate 
of the error bars on the observational data, both our result 
and the lognormal give a reasonable fit. 
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